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ABSTRACT 

The protoplanetary discs of T Tauri and Herbig Ae/Be stars have been studied by 
using geometric disc models to fit their spectral energy distribution (SED). The simula- 
tions provide means to reproduce the signatures of different circumstellar structures, 
which are related to different levels of infrared excess. Aiming to improve our pre- 
vious model that assumed a simple flat disc configuration, in the present work we 
adopt a reprocessing flared disc model that considers hydrostatic, radiative equilib- 
rium iDuUemond et al. 2001). We developed a method to optimise the parameters 
estimation based on genetic algorithms (GA). This paper is dedicated to describe the 
implementation of the new code, which has been applied for Herbig stars from the 
Pico dos Dias Survey catalogue, in order to illustrate the quality of the fitting for a 
variety of SED shapes. The star AB Aur was used as a test of the GA parameters esti- 
mation, demonstrating that the new code reproduces successfully a canonical example 
of the flared disc model. The GA method gives good quality of fittings, but the range 
of input parameters must be chosen with caution, since unrealistic disc parameters 
can be derived. The flared disc model is confirmed to fit the flattened SEDs typical 
from Herbig stars, however the embedded objects (increasing SED slope) and debris 
discs (steeper SED) are not well fitted with this configuration. Even considering the 
limitation of the derived parameters, the automatic process of SED fitting provides 
an interesting tool for the statistical analysis of the circumstellar luminosity of large 
samples of young stars. 
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1 INTRODUCTION 
1.1 Background 

The use of a model to study the disc geometry and 
the circumstellar emission from young stars is part of 
a detailed analysis that our group have developed for a 
large sample of T Tauri stars (TTs) and Herbig Ae/Be 
stars (HAeBes) detected by the Pico dos Dias Survey - 
PDS l|Gregorio-Hetem et al.i[l993 : [Torres et al]| 19951 : iTorred 
1 19981 ). In a general context, the interest is to discuss the 
geometry of protoplanetary discs in terms of disc evolu- 
tion. While the accretion processes occur for most of the 
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Among different circumstellar structures that can be 
used to explain the spectral ener gy distribution (SED) from 
young stars, in a previous work Gregorio-Hetem fc HetemI 
(2002) (hereafter GH02) used a flat disc model, which re- 
irradiates the energy absorbed from the star, leading to in- 
frared excess even in the absence of any accretion luminosity 
(lAdams et al ]|l987l : iKenvon fc Hartmannll 19871 : [ Strom et al.l 
Il993h . Assuming that flat disc models are valid when there is 
no accretion, GH02 adopted the simple model of passive disc 
to study weak- line TTs (WTTs). The occurrence of accre- 
tion is not expected for these objects, as indicated by their 
low Hq equivalent width (Whc)- They selected 27 P_DS stars 
showing Wua < 10 A, a usual criterion to distinguish WTT s 
from classical TTs (CTTs) (|Appenzeller fc Mundtl Il989l ). 
The individual contribution of the circumstellar components 
to the total emitted flux was compared with stellar charac- 
teristics in order to better classify the sample in agreement 
with an evolutionary sequence. GHQ2 concluded that almost 
half of the studied WTTs have considerable circumstellar lu- 
minosity (~ 30 percent of the total luminosity) confirming 
that the absence of accretion does not imply in the lack of 
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circum stellar material, as verified by lAndre fc Montmerlel 
l|l994l ) in a comparison between the near-infrared and Ha 
classifications of TTs. 

The contribution of the circumstellar matter to the 
SED slope is often used to recognize different categories 
of young stellar objects (YSOs) by following an observa- 
tion al classification based on the near-infrared spectral in- 
dex jLada fc W ilking"l984': 'Wilk ing et al.l 19891 : lAndre et all 
ll99S ). This classification suggests a scenario for the evolu- 
tion of YSOs, from Class to Class III, which is well estab- 
lished for TTs. 

Although the present work is focused on HAeBes, it is 
useful to describe the meaning of the different TTs classes 
in order to illustrate how the spectral classification relates 
to disc structure and evolutionary sequence. We summa- 
rize here the phases of pre-main sequence evolution, based 
upon the review by iFeigclson & Montmcrlc (1999) and fol- 
lowing the schematic illustration by ICarknej (|l998l ). The 
earlier evolutive phases are represented by Class (infalling 
protostar) and Class I (evolved protostar) . Both contain ob- 
jects that are surrounded by circumstellar shell and show 
increasing SED in the infrared range. Most of the Class II 
objects are CTTs, which are expected to have an accretion 
disc and to show flattened spectral index. The later phase is 
represented by Class III, in which are included WTTs and 
post-TTs. The disc is thin or non-existent in Class III ob- 
jects, since they have little or no infrared-excess and show 
steeper slope of the SED. 

For HAeBes the spectral classification and its relation 
with an evolutionary sequence are less advanced than the 
TTs classes. In Sect. ll.2l we present the HAeBes groups and 
a tentative scenario for their evolution, as proposed by other 
authors. 

1.2 The HAeBes groups 

The flat disc geometry was also adopted to model the 
protop lanetary disc s of H AeBes candidates detected by 
PDS. ISartori et ahl (|2003l ) (hereafter SGH03) classified 
these stars in different groups, based on the SED shapes 
and the fraction of circumstellar luminosity, expressed by 
fc = {Ltot — Li,)/Ltot . Adopting t he spectral index /3i = 
0.751og(Fi2/-Fv) - 1 (|Torres|[l99i '). the stars were classi- 
fied as follows: the higher infrared excesses (/c > 70%) are 
found in the /3i > group, which contains 44 per cent of 
the sample with increasing spectrum in the visible to in- 
frared range; the —1 < /3i < group has 50 per cent of the 
stars showing intermediate amounts of circumstellar emis- 
sion (/c = 10 to 70 %); in the third group (/3i < -1) there 
is 6 per cent of the stars having the lower values of infrared 
excess (/c < 10%). The groups were established following 
th e idea of an ev olutionary scenario for HAeBes proposed 
bv lMalfait et al.1 ([1998), similar to the TTs classes. In this 
scenario, values of j3\ > 0.7 correspond to embedded objects, 
like R CrA; and /3i < — 1 represents the stars showing more 
evolved discs, like Vega. The intermediary group has objects 
with SED similar to several known HAeBes. 

This classificati on is almost the s ame as the HAeBe 
groups presented bv lHillenbrand et al.l (|l992l ). based on the 
slope of the infrare d continuum. I t is al so comparable to the 
groups proposed bv lMeeus et al.l (|200ll ) that studied a sam- 
ple of 14 isolated Herbig stars. Their two major groups sug- 



gest different circumstellar geometries: the Group I sources 
are associated with an optically thin, flared region surround- 
ing an op tically th ick disc, while Group II sources have a flat 
disc ( Meeus et al.|[2002 ). One of the differences of this com- 
parison is related to the SGH03 group of embedded objects. 
Since they cannot be considered isolated Herbig stars they 
were not classified by Meeus et al. Another difference is the 
subdivision based upon the presence (a) or the absence (b) 
of the 10 /im silicate feature, which is not considered by 
SGH03. Table [1] summarizes (in the last columns) the com- 
parison of the groups above described. 

SGH03 estimated the circumstellar luminosity of the 
PDS HAeBes adopting the GH02 flat disc model. They ver- 
ified that the fitting-quality was gof > 0.4 (a bad data fit- 
ting) for more than 30 per cent of the objects. This suggests 
that different disc geometry should be more appropriate for 
an important part of the sample. 

1.3 Main objectives 

Wo decided t o adopt a fla red configuration, by following 
the DuUemond et al.l (|200ll ) (hereafter DDNOl) model for 
a passively irradiated circumstellar disc with an inner hole. 
Another improvement on our previous work is the method 
to optimise the parameters estimation based on genetic al- 
gorithms (GA). This automatic minimisation procedure 
was developed to increase the efficiency in the process of the 
SED fitting. 

The main goals of the present work are: (i) to describe 
the implementation of the GA method for the flared disc 
model; (ii) to discuss the derived disc parameters in a variety 
of SED shapes. Once the GA method has been conflrmed to 
be a valuable tool for SED fittings, the scope of the second 
goal of this paper is to discuss the value of different spectral 
classifications in determining disc structures; the relevance 
of flat versus flared geometries; and on degeneracy of the 
parameters derived from SED modelling. 

We organized the structure of the paper in the following 
sequence. In Sect.[2]we present the adaptation of the DDNOl 
model, for which a new code was developed. A test of the 
rewritten routines was made by applying them to AB Aur. 
Sect. [3] describes the GA method of parameters estimation 
and its implementation. Time efficiency and error estimative 
provided by the GA method are discussed in a comparison 
with our previous calculation procedure. In order to illus- 
trate the application of the GA method to fit different SED 
shapes, using the DDNOl model, we present in Sect. |4] the 
results obtained for some of the PDS stars. The sample was 
selected according different /3i ranges to represent different 
groups of HAeBe stars. We choose to study only a few stars 
in this first presentation of the GA method just to verify 
if the adopted disc geometry is suited to fit different SED 
shapes. In a forthcoming paper we intend to apply the GA 
method to evaluate the circumstellar structure of the PDS 
stars that probably have fiared discs. Considering the differ- 
ent levels of infrared excess shown by the selected sample, 
we discuss in Sect. [S] the circumstellar luminosity compared 
with the decreasing slope of the SED, which has been inter- 
preted as a sequence for the disc evolution. The quality of 
the SED fitting is used to compare different morphologies 
assumed in the flared and the flat models. For this purpose 
we describe in Appendix [Cl the GH02 flat disc model. Some 



Genetic algorithm to model protoplanetary discs 3 



outer l( 




+++++ 

+ + 



Figure 1. Geometric sclieme for tlie model adapted from DNNOl. 



inconsistencies of the derived parameters with respect to the 
observed ones are found for HD 141569, leading us to dis- 
cuss the degeneracy of the information obtained from SED 
modelling. Finally, in Sect.[6]the main conclusions are sum- 
marized. 



2 THE FLARED DISC MODEL 



IChiang fc GoldreichI (|l997l ) developed a flared disc model 
that assumes hydrostatic radiative equihbrium for passive 
dust discs of TTs. DDNOl proposed a modifled version of 
this model by including a puffed-up inner rim in the cen- 
tral disc region. This rim structure projects a shadow onto 
the disc and assumes that the disc's gaseous component is 
optically thin to steUar radiation. The dusty component is 
opticaUy thick, when the accretion rate through the disc is 
very low. 

To fit the SED we adopted the DDNOl model equations. 
For effects of calculation, the disc is divided into three parts: 
the inner rim, the shadowed region and the flared region 
(this has two layers: an illuminated hot-layer and an inner 
cold-layer). The first step is to establish a vertical boundary 
irradiated directly by the star, considering the effect due to 
shadowing from the rim, and the variations of scale height 
as a function of the radius. Figure [1] presents a reference for 
the geometry of the model. The FitCGPlus disc parameters 
are: radius, Rd, mass, Md, inclination, 9, density power 
law index, p, and inner rim temperature. Trim- The stellar 
parameters are: distance, d, mass, M^, luminosity, L*, and 
temperature, T*. 

We are aware that the treatment of the shadow by 
the inner rim in the DDNOl model is rather simple, and 
that full 3-D axis ymmetric radiative transfer models of 
iDuUemond fc Dom inik (2004) show that the shadow is in 
reality significantly less extreme. However, for the SED this 
does not appear to make too much of a difference. Another 
drawback of the DDNOl model compared to full radiative 
transfer models is the assumption that the inner wall ra- 
diation is a pure blackbody at the dust evaporation tem- 
perature: in 3-D radiative transfer calculations the inner 
wall is not a pure blackbody emitter and as a result pro- 
duces slightly weaker near-infrared emission (DuUemond, 
private communication). These differences are not expected 
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Figure 2. Comparison of the FiTCGPlus with the HABdisk re- 
sults for AB Aur. 



to strongly affect our results. The advantage of using the 
DDNOl model over full radiative transfer models is the speed 
of calculation of large numbers of models in a short time. 

In the present work we keep the original nomencla- 
ture and symbols described in Sect. 2 of DDNOl paper, 
which gives details on the physics of the model. C. P. DuUe- 
mond developed an interactive tool for applying the DDNOl 
model of protoplanetary discs. A public domain IDL code 
for their model, named FitCGPlus 2.1, is provided together 
with opacity tables for different grain composition, as well 
as photometric and spectral data of the star AB Aur, for 
testing. It is important to note that the IDL version does 
not reproduce the same level of details of the stellar spec- 
trum features (in the blue spectral region) as presented in 
Figure 8 of DDNOl. However, our results are not affected by 
the stellar spectrum, since our interest is mainly related to 
the SED fitting in the infrared region. 

Considering that our GA libraries are written using 
C++ language, and aiming to increase the speed of calcula- 
tion, we needed to rewrite the code for the DDNOl model. 
The code (called HABdisk) was built following a bottom- 
up strategy, and the result of each intermediary routine and 
function was compared to their correspondent in FitCGPlus. 

The comparison of the FiTCGPlus synthetic flux with 
HABdisk result gives a mean ratio of [(AFA)FiTCGPius — 
(AFA)HABdisk]/(AFA)HABdisk ~ 1-3%, that seems to be due 
to binary computations during iterations. Figure [5] illus- 
trates these small differences by showing the deviations for 
each considered wavelength. 

The HABdisk code was applied to AB Aur (a DDNOl 
canon ical model) using the parameters from lDominik et al.l 
(|2003l ): = 2.5 Mq, L* = 47 Lq, T* = 9750 K, Md = 0.1 
Mq, Rd ^ 400 AU, 6^ = 65°, Trim = 1500 K, and D = 144 
pc. Figure [3] presents a comparison of the SEDs calculated 
with the codes FitCGPlus and HABdisk. 



3 THE GA METHOD 

This section explains how the GA method was implemented 
to the DDNOl model (more details can be found in the Ap- 
pendix [XJ. 

Considering that the GA nomenclature can be mis- 
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Figure 3. AB Aur synthetic SEDs obtained witli two codes. The continuous thin Une represents the star emission; dashed Une represents 
the rim emission; dot-dashed line for the disc cold-layer; and dotted-line for the disc hot-layer, (left) FitCGPlus, the IDL version of 
DDNOl model; (right) HABdisk, the C+-f version of DDNOl model. The plot is given in log [aF;^ (Wm~^)] versus log [A (fim)]. 



placed in the astrophysical context, we first present the 
translation from one field to anotfier. A parameter (e.g. disc 
radius) , corresponds to tfie concept of a "gene" , and a change 
in a parameter is a "mutation" . A parameters set tliat yields 
to a possible solution corresponds to a "chromosome". An 
"individual" is a solution, that is composed by one param- 
eters set and two additional GA control variables. One of 
these variables is i which means the "adaptation" level. 
The other control variable is <E>, the genetic operator de- 
scribed in Sect. I3.IT2] The term "generation" means "all the 
individuals" (or all the solutions) present in a given itera- 
tion. 

3.1 Implementation 

The HABdisk code uses the same parameters of FitCGPlus, 
which are: M*, Rd, 6, Md and p. On the other hand, we 
decided to allow Trim to be variable. The stellar parameters 
d, and T*, are adopted from observations (see Sect. 14.2)) . 

Essentially, the GA method presented herein imple- 
ments a minimisation of the SED fitting provided by 
the DDNOl model. There are three main advantages in us- 
ing GA for this task: (i) the GA method potentially browses 
the whole permitted parameter space, better avoiding the 
"traps" of local minima; (ii) the method is not affected by 
changes in the model; (iii) the GA implementation does not 
need to compute the derivatives of (hke dx^/dRo for 
example) required by the usual methods. This fact simpli- 
fies the code and minimises computer errors due to gradient 
calculations. 

3.1.1 Data structures and control variables 

The main structures used to manipulate the data are linked 
lists containing the solutions (parameters set, the adaptation 
level, Xiy a-nd the genetic operator, expressed by 

S,^ [{RD^,e,,MD^,P^,T) ,{xl^^)] W 

where 5"; mean s the zth s o lution , and Ti is the ith Trim. 

FoUowing iGoldberd (|l989l ) the code starts with the 
construction of the first generation, where all parameters 



are randomly chosen within an allowed range (for example, 
50 Rd ^ 1000 AU). In the present work, the number of 
different parameter sets in the first generation is assumed 
to be 100. In the next step the evaluation function runs the 
DDNOl model for each solution, and compares the synthetic 
SED with the observed data to fi nd the x^, using the expres- 
sion given bv lPress et all l| 19951 ): 

N 

X? = ^^(i^.-0..)' (2) 
i 

where A'^ is the number of observed data points, Fj is the 
observed flux at wavelength \j and (j)ij is the calculated 
flux for the solution 5*;. The smallest x^ corresponds to the 
goodness-of-fitting, or gof. 



3.1.2 The genetic operators 

As the evaluation function has been applied to all solutions, 
the judgement procedure sorts the list by increasing x^ , and 
sets to "l>i one of the genetic operators: copy, crossover (also 
called recombination), mutation or termination. Each $ is 
attributed to a f raction of th e generation follow i ng the val- 
ues suggested bv lKozal (|l994l ). lBentlev fc Cornj l|2002l ) and 
references therein. 

The next generation is constructed by applying the cor- 
respondent genetic operators. The copy uses an elitist se- 
lection, since the solutions with the smallest Xa £^re copied 
to the next generation. The crossover operator randomly 
mixes parameters from two different solutions. For example: 

Sa = [RDa,Oa,MDa,Pa,Ta\ aud Sh = [Rob , db , M Db , Pb ,Tl,] 

can give Sc = [RDb,Oa, MDa,Pb,Ta]. Mutation reproduces 
the same original solution, except for one of the parameters 
that is randomly changed. 

The method continues to build new generations until 
the end condition is reached, as illustrated by the schematic 
view shown in Figure lAll The distribution of x^ for ^ ^6" 
quence of 50 generations in the SED fitting of AB Aur is 
presented in Figure [4] 
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Figure 4. The distribution of for a sequence of generations. 
The circles are used to represent the best solution in each gen- 
eration, and the triangles represent the worst ones, the 5% best 
solutions are found below the level indicated by the crosses. 
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Figure 6. gof contour levels in parameter space projected in 
[Af^i, Ro] plane evaluated for AB Aur. Levels are = 0.082 
(continuous line), = 0.15 (dashed) and = 0.21 (dot- 
dashed). 
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Figure 5. Histograms of x^ values, comparing two different min- 
imisation methods used to fit the SED of AB Aur. The Newton- 
Raphson method results are filled with horizontal hatching, while 
the GA method results are filled with 45° hatching. 

3.2 Tests and error estimative 

3.2.1 Comparing to previous method 

Figure [5] shows a comparison between the GA method and 
the classical Newton-Raphson method adopted in our pre- 
vious work (GH02), both used in the SED fitting of AB 
Aur. The test is limited to 256 runs, each one starting from 
random values for the free parameters. The number of itera- 
tions is fixed, aiming to obtain a time independent compar- 
ison. The GA method provides a single concentration that 
approaches the global minimum around = 0.046. The 
spurious peaks produced by the classical method are due to 
the several local minima that usually result from methods 
based on gradients. We also evaluated the time efficiency 
and verified that the GA method is about 16% faster. 

3.2.2 Confidence levels 

An analysis of the behaviour as a function of the pa- 
rameters variation is used to esti mate error bars, de termin- 
ing the confidence levels given bv lPress et al.l (|l995l ). as de- 
scribed in AppendixlB] The error bars estimated for la con- 
fidence level and respective disc parameters for AB Aur arc 



Md = 0.1 ± 0.004 Mq, Rd = 400 ± 44 AU, 5 = 65 ± 3°, 
and Trim ~ 1500 ± 26 K. These results are in agreement 
with another error bars estimation provided by the surface 
contour levels. 

Figureinipresents the contour levels of the gof{MD,RD) 
surface, which was calculated for a set of 400 random pairs 
of disc mass and r adius around the para meters for AB Aur 
model taken from [Pominik et al.l (|2003l ). resulting a mean 
minimum gof ~ 0.046. This means that the GA method 
converged to a narrow range around the same parameter set 
obtained by them. The three contours represent the confi- 
dence levels for 68%, 90% and 99%. 



4 TESTS WITH PDS STARS 

4.1 Sample Selection and observational data 

Motivated by the efficiency of the GA method, in a further 
work we intend to model the protoplanetary dis cs of the 108 
HAeB e detected by PDS, which are listed by IVieira et al.l 
l|2003h . The worth in evaluating the circumstellar luminosi- 
ties of a large sample of objects is to compare them with 
spectral and polarimetric characteristics, providing a statis- 
tical analysis to be presented by our group in a forthcoming 
paper. The main goal of the present work is to apply the GA 
method for a few objects only, in order to verify the qual- 
ity of the fitting for objects showing different SED shapes 
and different levels of infrared excess. We decided to choose 
four stars to represent the groups described in Sect. 11.21 
selected by the slope of their S ED in the near-infrared. Ac- 
cording Ivan den Ancker et al.l l|200ll ). the infrared excess in 
HBe stars is due to a spherical dusty envelope, while HAe 
stars are better represented by a thick edge ffared disc. This 
suggestion motivated us to select A type or late-B type stars 
from the PDS sample to illustrate the application of the SED 
fitting procedure described in Sect.|3l Table [1] presents the 
main characteristics of the stars considered representatives 
of different ranges of f3i discussed by SGH03. 

There is a lack of information about IRAS 07394-1953 
and IRAS 6475-0735, being only previously studied by th e 
Pi?5 group (|Gregorio-Hetem et al. I ll992l : IVieira et aLll2003l ) . 
BD-14 1319 has been included in a polarimetric study by 
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Table 1. List of P_D5 stars selected as representatives of different /3i groups (SGH03 ). Spectra l types , estimated visu al extinctions and 
stellar parameters are also given. A comparison with HAeBe groups suggested by Me eus et alj 1I2OOII) and [Hillenbran d et al.l lll992h are 
presented respectively in the two last columns. IRAS 07394-1953 (PDS257) is not classified in the groups of Meeus et al., since it cannot 
be considered an isolated Hcrbig star. 



PDS 


Name 


S.T. 


Av 
(mag) 


(K) 


(Lq) 


d 
(pc) 


/3i 


SGH03 


Meeus 
group 


Hillenbrand 
group 


398 


HD 141569 


B9.5Ae 


0.25 


10000 


20 


106 


-1.74 


/3i < -1 


II 




022 


BD-14 1319 


AO 


0.37 


9760 


78 


800 


-0.42 


-1< /3i < 


I 


I 


130 


IRAS 06475-0735 


A 


1.9 


10140 


35 


995 


0.23 


< /3i < 0.7 


I 


II 


257 


IRAS 07394-1953 


AO 


2.7 


10000 


150 


4000 


0.74 


/3i > 0.7 




II 



lYudin fc Evand (|l998l ) but no polarization was detected for 
this star. Complementary p hotometry and low-re s olutio n 
spectroscopy are reported bvlMiro shnichenko et al.l (Il999l ). 

On the other hand, HD 141569 is a very well studied 
object, showing a circumstellar structure that has been con- 
sidered similar to the disc of /3 Pictoris. The B9.5Ve primary 
and two T Tauri stars form a triple syste m, separated by 9 
arcsec , with age of 5 Myr, as estimated bv lWeinberger et al.l 
(|2000l ). iMen'n et all (|2004l ) estimate L = 22 Lg for the pri- 
mary. Based on the study o f circumstellar accretion of a 
sample of PDS HAeBe stars, iGuimaraes et all l|2006l ) have 
determined stellar parameters, metallicity and rotational ve- 
locity. No evidence of accretion episodes was found for HD 
141569, probably related to an evolved status of the disc. 

A NIC MOS image of the circu mstellar structure was 
obtained by [Weinberger et al.l (|l999l ) showing a disc size of 
almost 5 arcsec (500 AU in projection). Observations of the 
disc in scattered light at near-infrared and optical ranges re- 
vealed a complex dust distribution showing a large inner hole 
out to 150 A U from the central star and two rings at 185 AU 
and 325 AU ('Augcrcau et al.lll9 99': 'Weinberger et al."l99gl; 
Mouillet et al. 2001; Boccaletti et al. 2003_; Clamp in et al.i 
20031 ). 

The r emnant mass of H2 w as estimated as 20-460 
Mgarth by IZuckerman et al.l (Il995l ) based o n CO measure- 
ments around HD 141569. More recently Ijonkheid et al.l 
l|2006l) have studied the chemistry and gas temperatures in 
the disc using a radiative transfer code to fit the observed 
CO emission lines. In order to reproduce the data, they 
adopted a power-law gas distribution of 80 Msarth start- 
ing at 80 AU from the central star. 

The optical photometry was obtained from PDS and 
far-infrared from IRAS Point Source Catalogue. An up- 
per limit millimetric photometry at 1.1 mm, obtained by 
ISvlvester et~aLl l|l996l ) is available only for HD 141569, which 
gives a constraint for the cold layer disc component. Near- 
infrared data for IRAS 06475-0735, IRAS 07394-1953, and 
BD- 14 1319 are fror t i 2M ASS catalogue, and HD 141569 
from ISvlvester et"al] (|l996l ) . 



ISchultz fc Wiernerl (Il975|). The intrin s ic co lours and bolo- 
metric corrections from iBessel et al] (|l998l ) were adopted 
in th e estimation of Ay and absolute bolometric magni- 
tude. IVieira et al.l (|2003i ) have listed the association of the 
objects with the star forming regions CMaRl (1000 pc, 
IRAS 06475-0735), NGC2149 (800 pc, BD-14 1319), and 
G236.4-f-l.49 (4300 pc, IRAS 07394-1953). For these stars 
we have adopted the distances of the respective clouds and 
assumed an error of 10%. Hipparcos parallax has been ob- 
tained for HD 141569 giving a better determination for dis- 
tance (99 ± 10 pc) in this case. 

Temperatures were estimated from the T^ff- Spectra l 
Type relations given by lde Jager fc NieuwenhuiizenI l|l987h . 
assuming a deviation of T* ± 200 K. 

Table [T] gives the distance, luminosity and temperature 
of the stars that are fixed input parameters. The free param- 
eters and respective adopted limits are: stellar mass (2 to 20 
Mq); disc mass (0.001 to 0.5 Mq); disc radius (100 to 1000 
AU); and inclination angle (0° to 70°). The only exception 
is HD 141569, for which we have fixed the disc size as 500 
AU. 

Two of the model parameters ar e not fitted by u s in th e 
same way as made by DDNOl and iDominik et al] l|2003h . 
The slope of the power-law E(7?) = EqT?^, which describes 
the surface density, is assumed not to be a completely free 
parameter, since only the discrete values p= -2, -1.5 or -1 
were used. On the other hand, the disc inner temperature is 
not fixed at 1500 K. Instead of fitting the inner disc radius, 
we adopted Trim as a free parameter, ranging from 300 to 
2000 K, in order to have different si zes of the inne r disc h ole. 

The dust opacity given by iDraine fc Led (|l984l ) is 
adopted, assuming a 100% silicate grain composition. It 
must be kept in mind that different grain compositions cause 
significant degenerations in the model solutions (see Sect. 
I5.1|l . We are not attempting to fit in detail the spectral fea- 
tures from the near- to the far-infrared, considering the lack 
of detailed observations, like ISO data, for most of the PDS 
stars. For this reason, we decided to adopt a standard grain 
composition for the whole sample. 



4.2 Input parameters 



4.3 Results 



The optical and near-infrared data were corrected for in- 
terstellar extinction b y using the normal rel ation Ay = 



3.1E{B - V) given b ylSavage fc MathisI ( 19791 ) and extinc- 
tion law A\/Av from lCardelli et al.l ( 198E ). Considering the 



UV excess presented by Herbig stars, the {B — V) excess 
was estimated by adopting E{B — V) — E(V — /)/1.6 from 



Table [5] presents the parameters derived from SED fitting 
by the HABdisk code, and the fraction of circumstellar lu- 
minosity (/c). A very good quality of fitting was achieved 
for BD-14 1319 and IRAS 06475-0735, as shown in FigureEl 
The same cannot be said for the other two objects displayed 
in Figure [8] 
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Table 2. Derived model parameters for the sample, where fc gives the per cent amount of circumstellar luminosity. The best fitting of 
HD 141569 is provided by parameters set (b). 



PDS 


Name 




(Mq) 


Rd (AU) 


Md (Mq) 


i7n ( ) 


9 n 


P 


gof 


fc (%) 


/3i 


398 


HD 141569 


(a) 


2.4 


500 


0.002 


303 


0.6 


-2.0 


0.021 


1.4 


-1.74 






(b) 


2.4 


13 


0.06 


1085 


0.6 


-2.0 


0.006 


1.4 


-1.74 


022 


BD-14 1319 




2.8 


690 


0.003 


380 


40 


-1.0 


0.006 


21 


-0.42 


130 


IRAS 06475-0735 




2.0 


309 


0.20 


1705 


53 


-1.5 


0.016 


32 


0.23 


257 


IRAS 07394-1953 




2.0 


859 


0.64 


1838 


47 


-2.0 


0.098 


> 35 


0.74 




Considering the low quality of the fitting obtained for 
HD 141569 when fixing Rd = 500 AU, we decided to test a 
second HABdisk run for this star, and to compare different 
sets of parameters. In this case, an excellent SED fitting was 
obtained (see Figure but the derived disc parameters are 
unrealistic, a s shown in T able [2] (see also Sect. 15.3")) . 

Domini k et al.1 (120031 ) have discussed the variety of the 
spectra of is o lated HAeBes, members of groups I and II of 
iMeeus et all l|200ll ), in the context of the DDNOl model. 
They found good SED fittings for both groups, especially in 
the case of Ila objects. However, some of the fitting resulted 
in compact disc radii that are inconsistent with the typical 
sizes measured around young stars. They concluded that 
while the DDNOl model is capable to fit the SED of Group 
Ha objects, the derived model parameters seem difficult to 
understand. 



5 DISCUSSION 

5.1 The quality of the fittings 

In this section we first discuss the results of goodness-of- 
fitting and the adequateness of the DDNOl model for the 
studied objects. It is confirmed that the flared geometry is 
well suited for the stars showing intermediary values of Pi 
index, while the same cannot be said for the objects with 
extreme SED slope (/3i < -1 or /3i > 0.7). 

Despite the good fittings, only IRAS 06475-0735, which 
shows a fiattened SED similar to AB Aur, had more realis- 
tic derived parameters. For BD-14 1319, an unrealistic Trim 
was required to fit its double peaked SED. A large central 
gap yielding t o a coo ler inner disc temperature, as found by 
lAlencar etal] (|2003l ) in the SED fitting of AK Sco for exam- 
ple, could be the reason. However, the derived Trim is much 
lower than the usually adopted dust evaporation tempera- 
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2 



Figure 9. Best SED fitting obtained for HD 141569, which gives 
Rd = 13 AU (see text). 



ture. The absence of a pufTed-up rim seems to be a better 
explanation, since this star does not present a bump in the 
near-infrared. 

It can be noted from Figure [8] that the object showing 
the largest infrared emission, IRAS 07394-1953, could not 
be fitted by the adopted model. A possible reason is con- 
tamination by nebulosity or cirrus as indicated by the high 
excess at 60 nm and 100 /im verified on the IRAS-ISIS im- 
ages, which respectively show flux densities Fqq > 5 MJy/sr 
and -Fioo > 20 MJy/sr nearby the location of IRAS 07394- 
1953. For the other fields, lower levels are found {Fao ^ 3 
MJy/sr and Fioo ^ 10 MJy/sr), leading us to suggest that 
IRAS 07394-1953 is the most embedded object of the sam- 
ple. Probably in this case, a circumstellar shell component is 
needed to reproduce the far-infrared excess, as discussed in 
the comparison of different models (see below). Only a lower 
limit of fc > 35% could be estimated for this star since the 
infrared data were not fitted. On the other extreme, HD 
141569 does not have a good fitting if the observed disc size 
is adopted (see Sect. I5.3"|l . 

In order to verify the question of whether a flat disc 
model would provide good SED flttings for the HAeBes stud- 
ied in the present work, we also present previous results ob- 
tained with the GH02 model, which is described in Appendix 
[Cl Table ICll gives the parameters derived for this sample, 
showing that the quality of the fitting (gof) is worst than 
that obtained with HABdisk, except for IRAS 07394-1953. 
In this case, the spherical component included in the flat disc 
model could provide a better fitting of the far-infrared data. 
Independently of the adopted circumstellar geometries, the 
poor gojs achieved with the flat disc model seems to be re- 
lated to the minimisation procedure, as discussed in Sect. 
13. 2. H and shown in Figure [5] 

Since a flat disc would be expected to explain the HD 
141569 circumstellar structure (/3i < — 1, Group II of Meeus 
et al. 2002) we compare here the flat disc parameters with 
the HABdisk results, obtained for this star. The disc size 
and the temperature resulted in reasonable values using the 
flat disc model, but the system is seen almost edge-on, con- 
tradicting the observations. Probably, this is due to the high 
disc luminosity produced by the optically thick disc assumed 
in the GH02 model, as expressed in Eq. IC3I 

The differences and similarities resulting from both 



Table 3. Disc parameters obtained with the SED fitting of HD 
141569 by adopting different grain compositons. 





Rd (AU) 


Md (Mo) 


im { ) 




carbon 


21 


0.012 


1050 


0.1 


forsterite 


24 


0.012 


748 


0.04 


H2O ice 


65 


0.010 


576 


0.01 




og[A (urn)] 

Figure 10. HD 141569 ISO photometry (filled circles) from SS 
detector (2.5-4.9 ^m) and SL detector (5.8-11.6 fim), and SWS 
spectroscopy. Several curves are used to represent the SED fittings 
obtained with different chemical compositions. 



models demonstrate that degenerated information can also 
be derived from SED modelling. In this case, there is no ar- 
gument to distinguish between HABdisk and the GH02 flat 
disc model, when applied to HD 141569. 

Another consideration regarding the problem of degen- 
eracy is to verify how much grain composition affects the 
synthetic spectrum. For this study we have inspected the 
features shown by ISO photometry and spectroscopy avail- 
able for HD 141569 (public archival data). The details in 
the spectral range 2 to 50 /im are presented in Figure [10] 
that compares the ISO data with the SED calculated with 
the GA method by using different grain compositions. Table 
[3] presents the obtained disc parameters and gofs. Among 
the three grain types, carbon and forsterite gave the worst 
gofs, when compared to the results for HD 141569 in Table 
[2I On the other hand, H2O ice provided a good quality of 
fitting, but the mid- and far-infrared data were not well re- 
produced. The derived parameters are in the range defined 
by the results (a) and (b) presented in Table [2] indicating 
that the conclusions about HD 141569 are not biased by the 
choice of using astronomical silicates in the modelling. 



5.2 The correlation of fc and SED slope 

The four objects of the sample show a correlation of cir- 
cumstellar luminosity with spectral index as can be 
noted in the two last columns of Table (2] This correlation 
defines a distribution of the objects, where IRAS 07394- 
1953 (/3i = -hO.74, fc > 35%) is in the beginning and HD 
141569 (/3i = -1.74, fc = 1.4%) is in the end of the se- 
quence. The diminishing of fc with the SED slope is in agree- 
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Figure 11. H-R diagram sliowing tlie position of the stars com- 
pared witli the ZAMS (cont inuous line) and 10, 5 and 1 Myr 
isochrones (dotted lines) from lSiess et alj l|2000l ) . Only the largest 
error bars are presented. 



ment with th e HAe Be evolutionary sequence suggested by 
iMalfait et all (|l998l ). 

In a comparison with the circumste llar geometries of th e 
two major HAeBes groups proposed bv lMeeus et al.ll|2002h . 
the low fc of HD 141549 is consistent with their suggestion of 
a fiat disc structure for Group II, while the intermediary fc 
values for BD-14 1319 and IRAS 06475-0735 are consistent 
with the optically thick disc surrounded by a flared thin re- 
gion, as proposed for their Group I. On the extreme is IRAS 
07394-1953, showing the highest fc of the sample, which pre- 
cedes the two groups of Meeus et al. in the sequence, since 
it is not an isolated HAeBe. However, there are no clear 
differences among the derived model parameters for these 
stars that were chosen to represent different groups, or SED 
shapes. Based only on the SED fitting it is difficult to dis- 
tinguish the objects in an evolutionary sequence, mainly for 
those showing — 1 < /3i < 0.7. 

Indeed, the comparison of the stellar parameters on the 
H-R diagram shows no evidence of differences on their evolu- 
tionary status. In Figure [TT]the bolometric luminosities and 
effective temperatures of t he four sta r s are compared with 
the isochrones provided bv lSiess et al] l|2000 !). The Zero Age 
Main Sequence (ZAMS) and isochrones for 10, 5 and 1 Myr 
are displayed. Even considering the well-determined age of 
HD 141569 in the literature, it is not possible to infer the 
ages of the other objects in Figure [TT] due to the imprecision 
on distance or temperature determination, as indicated by 
the illustrative error bars. In spite of the age uncertainties, 
it can be noted that high values of /3i correspond to higher 
positions above the ZAMS, but a large number of objects is 
needed to verify this trend. 

In fact, the only hypothesis for the studied sample is the 
suggestion regarding the extremities of a possible sequence. 
Pi < —1.0 objects showing very little circumstellar emission 
(/c < 10%), as the example of HD 141569, are the best can- 
didates to have more evolved discs, while IRAS 07394-1953 
probably is the opposite case. Any other inferences about 
the evolutionary status of objects require a demonstration 
for a large sample that we intend to develop in a further 
work. 



5.3 HD 141569 

For HD 141569 the SED could not be well fitted with the 
disc radius measured from NICMOS image, indicating that 
the flared disc model fails to give the best fltting mainly in 
the mid-infrared region around ~ 25 fim. This is not due to 
a failure of the GA method to flnd the real best solution. 
It is verified that a different parameters set can provide an 
adequate SED fitting, but resulting in disc radius Rd ~ 13 
AU. While this seems to contradict the near-infrared obser- 
vations, which revealed a disc size larger than 300 AU, it 
must to be noted that only sub-/xm sized grains are seen in 
scattered Hght at large distances. The outer parts of the disc 
do not comprise the bulk of the disc mass, which is mostly 
concentrated closer to the star. 

Probably, the evolved disc (or thin disc geometry) of HD 
141569 is better explained by a flat structure instead. Our 
previous results obtained using the fiat disc model are con- 
sistent with the observed size (see Appendix [Cj. However, 
that simple model also fails to reproduce the complicated 
structure of the inner disc of HD 141569 that seems not to 
be only related to differences on circumstellar geometry or 
degeneracy of the models. We trace below some informa- 
tion about grain sizes and chemical composition, probably 

related to evolutionary aspects. 

The pr esence of large grains jBoccaletti et aL 



1 he pr esence oi large grams llpoccaletti et al. 
and PAHs jlSvlvester et al.1 Il996l : [Weinberger et al 



iGeers et al.]|2006l) have been repor ted for HD 141569. whil e 
silicate emission is not present (jWeinberger et al] |2000| ) . 
This can be verified in the spectral features presented in 
Figure [TOl The 750 data show an almost flat SED near 10 
^m, while a low intensity PAH emission is found at 8.6 /im. 

The absenc e of th e 10 fim feature has been discussed 
by iMeeus et al.1 l|2002l) as a possible lack of small silicate 
grains, or a temperature effect that leads to a lack of small, 
hot silicates. According DDNOl, the shadowing caused by 
the rim can suppress t he strength of the silicate emission. 
[Hernandez et all (|2007l ) suggest that grain growth and/or 
settling, and transitions objects, which are stars w ith in- 
ner gaps in their discs (|Sicilia-Aguilar et al.l |2006|) , are a 
possible combination to explain evolved discs. The correla- 
tion bet ween the mid-infra red and the sub millimetre slope, 
found bv lAcke et al.l (|2004l ) , indicates that as the grains grow 
the disc structure evolves from flared to geometrically flat 
(|Dullemondl2002l ), which seems to be the case of HD 141569. 

Another implication on the radial structure of this 
disc can be also related to dynamical processes. Scattered 
light images of the op ticall y thin dust disc were u s ed by 
lAugereau et al.l (Il999l ) and lAugereau fc Papaloizoul (120041 ') 
in the investigation of the properties and the dynamics of 
the HD 141569 triple system. Their suggestions on different 
dust populations required to fit the full SED, as well as the 
processes of disc truncation, indicate a scenario too complex 
to be reproduced by the two models that we have studied. 



6 SUMMARY AND CONCLUSIONS 

The goal of this paper is to describe a method to fit the SED 
of pre-main sequence stars, using the code called HABdisk, 
which was adapted from the flared disc model proposed by 
DDNOl. The procedure uses the GA method in order to find 
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the best parameters set, based on minimum criterion. 
The tests apphed to AB Aur confirm that the GA method 
provides a good estimation of the circumstellar parameters, 
reproducing a canonical DDNOl model. 

Aiming to provide an illustrative application of HAB- 
disk, we selected four PDS stars that have similar mass and 
effective temperatures, but show different amounts of in- 
frared excess, as indicated by their spectral indices. 

6.1 The value of the different spectral 
classification in determining disc structures 

The estimated circumstellar luminosities are in agreement 
with the expected correlation of SED shape with /c, which 
vanishes as the spectral index in the infrared decreases. The 
fc behaviour in the sample is consis tent with the Groups I 
and II defined bv lMeeus et al.l l|2002l ). as summarized below. 

The intermediary /c values estimated for BD-14 1319 
and IRAS 06475-0735 are in agreement with the optically 
thick disc surrounded by a flared thin region proposed for 
Group I. 

HD 141549 shows the lowest value of /c that is consis- 
tent with the flat disc structure suggested for Group II. 

Among the four stars, IRAS 07394-1953 is the only one 
associated with nebulosity, as verified in the IRAS-ISIS im- 
ages (see Sect. 15. 1[) . and probably for this reason it has a 
rising SED at the far-infrared that could not be fitted by 
HABdisk. 

6.2 Degeneracy of information from SED 
modelling 

Despite the good quality of fittings, our results also indicate 
that some degenerated information from SED modelling can 
occur, as verified in the comparison of different models. 

First, the fiared disc model was not able to provide a 
good SED fitting when the observed disc size was considered 
for HD 141569. A bad performance of the GA method is not 
the case, since a good fitting was achieved in a second run 
of the code, but leading to unrealistic disc parameters. 

Considering that a flat disc would better explain the cir- 
cumstellar structure of HD 141569, we reanalyse the results 
from GH02 model, which gave reasonable disc size and tem- 
perature. However, other unrealistic parameters were also 
derived, probably due to the adopted optically thick regime 
that gives high disc luminosity (see Table [CT|) . This is not 
compatible with the optically thin disc that is observed in 
scattered light at large distances from the central star. The 
disc parameters determination from the SED fitting remains 
unsolved since the radial structure of this disc is much more 
complicated than it is possible to describe using the adopted 
simple models. 

6.3 Flared vs. flat geometry 

On the view of the good results found for IRAS 06475-0735 
and BD-14 1319, we conclude that PDS stars within the 
-1 < /3i < 0.7 range have SEDs better fitted by HABdisk 
instead of the fiat disc model. 

For objects with /3i > 0.7 another component seems to 
be required in the circumstellar structure, such as a shell for 



example, to better reproduce the observed excess at longer 
wavelengths. In this case, our previous flat disc model, which 
includes a thin dust envelope in the system, provided better 
results than HABdisk. 

Even considering the degeneracy in the case of HD 
141569, its low fc value lead us to expect that objects show- 
ing /3i < — 1 seem to be better explained by flat disc geom- 
etry models. 

6.4 Conclusions 

On the light of the above remarks, we conclude that good es- 
timations of the circumstellar luminosities may be provided 
for a large sample of objects, and a statistical analysis can 
be done by comparing fc with /9i in order to verify a possible 
sequence of disc disappearance. 

We are aware that no inferences about disc evolution 
can be attended from the derived parameters. Considering 
the lack of detailed observations about disc parameters or 
spectral features related to grain composition, a parameters 
set must be carefully evaluated, even when derived from 
good SED fittings. Since these detailed observations are not 
available for the large majority of the PDS sample, the disc 
parameters must be constrained to realistic ranges. 

The difficulty to infer disc parameters for large samples 
arises from the interpretation of details on infrared spectral 
features and/or high angular resolution data, requiring in- 
dividual analysis of the physics, the geometry and chemical 
composition. This task shall be accomplished with the use of 
several promising high-sensitive observational facilities that 
are under operation or will be in a near future. In particu- 
lar, we mention here the near-infrared coronagraphic imager 
(NICI) that is being used with Gemini telescope to search 
for exoplanets. One of the sub-products of this campaign cer- 
tainly shall be the detailed imagery of protoplanetary discs. 
Using the GA method, we seek to select good targets among 
the PDS stars that are in the distance range to be imaged 
with NICI. 
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APPENDIX A: GENETIC ALGORITHMS 

GAs are a family of computer models inspired on natural 
evolution. Their basis is to assume a potential solution for 
a specific problem, viewed as chromosome like structure, on 
which are applied genetic operators (mutation, crossover, 
adaptation and evolution). The use of this technique simpli- 
fies the formulation and solution of optimisation problems, 
and parallel simultaneous procedure approach is implicit in 
the method, providing evaluation of the viability of a pa- 
rameter set as possible sol ution for complex problems ijKozal 
1 19921 . 1 19941 : lHoUand|[l993 ) . 
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Figure Al. Schematic view of the main steps of a generic GA. 



The implementation of a GA starts with a random gen- 
eration of sets of chromosomes [Si, S2, ■ ■ ■ , Sn], which are 
evaluated and associated to an adaptation probability, x^y 
obtained by the evaluation function, which is in our case 
the DDNOl model. The Xi values express how each individ- 
ual are a dapted, or how each so lution are close to the best 
solution l|Bentlev fc Cornell2002D . 

Then, the judgement function determines the genetic 
operator to be applied to a solution, and its values can be 
copy: the individual remains the same in the next genera- 
tion, crossover: the individual is elected to change a num- 
ber of genes (parameters) with another individual, creating 
a new one, mutation: one of its genes shall be randomly 
changed, or termination: none of the genes may continue on 
next generations. The chosen action is expressed by the "l>i 
variable, associated to each individual. The further step is 
to evolve the current generation (fc) to the next (fc + 1), what 
is done through a procedure F that considers the solutions 
and the genetic operators designed by <l>i. Formally 



[Si, S2, .., S„] 



k + l 



r[(5i,<&i),(S2,$2),..,(S„,<&„)], (Al) 



As soon as a new generation is ready, the evaluation 
function is reapplied, and the algorithm repeats the de- 
scribed actions until an end of loop condition is reached. 
The end condition can be based on the number of iterations 
or quality (a low level for the Xi values) (Figure TAIII . 



APPENDIX B: ERROR BARS ESTIMATIVE 

Once an acceptable minimum has been found, one computes 
the inverse of the Hessian matrix 



[C] ^ [a]-' 

that has its components given by 
dy{\k) dy{Xk) 



N 

E 



dui daj 



(Bl) 



(B2) 



where dy{Xk)/dai means the partial derivative of the SED 
with respect to parameter at A = Afc, and is the number 
of observed data points. The components of the main diag- 
onal of C can be used to estim ate the error bars on each 
parameter by l|Press et al.lllQQSI '): 
^1/2 

^ ^ (B3) 



Table CI. The SED fittin g parameters obtained by applying the 
flat disc model ijGregorio-H ctcm & Hctcm 2002) for the sample. 
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APPENDIX C: FLAT DISC MODEL: REVIEW 

GH02 used a flat disc model assuming that the star and the 
disc are surrounded by a thin dust shell. The emission from 
the star and from the disc are attenuated by the opacity of 
the envelope. The flux is calculated by assuming blackbody 
emission and different temperature laws for the disc and the 
envelope. 

The model considers the following free parameters: ra- 
dius Ra of the star, radius Rd and inclination 9 of the disc, 
radius Re and total optical depth r of the envelope. The 
stellar temperature Ts is a fixed parameter, since it can be 
derived from the spectral type. The absorption and the emis- 
sion of radiation depend on t, where the extinction law is 
normalized to r(l/im) = 1. The stellar flux is expressed by 



Ss{\) = ^ttB (A,r.)e-^(^) 



(CI) 



A passive, optically thick and geometrically thin disc is 
assumed, for which the radia tion is evaluated by adopting 
the temperature law given bv I Adams Sz Shul (|l986l ): 



-3/4 



(C2) 



The inner radius of the disc is defined by Rh = 
Rs{Tg/Ts)~*^^, with adopted Tg =1500 K as the destruc- 
tion temperature of the grains, by consider ing the silicate 
sublimation temperature (jChiang et al.ll200ll ). The disc con- 
tribution to the total emitted flux is given by 



2r5r 



ttB [A,T(r)]e"^<^' cos 6'dr 



(C3) 

the radiative transfer model by 



Sd (A) = 

Following 

lEpchtein etall l|l990t ). the emission from the dust en- 
velope is calculated using the temperature law 

2r 



Te{r) = T, 



(f) 



(C4) 



The optical depth depends on the density distribution: 
T(r) oc p(r), where p(r) oc r~^^^. The flux emitted by the 
envelope is obtained by 



Se{\) 



B[\,Te{r)]e 



-T(r,\) 



-Te(r,A)N 



dr(C5) 



where is the optical depth calculated from r to J?e- 

The total radiation emitted as a function of wavelength 
is expressed by ^(A) = ^^A) + Sd{\) + Se{\). 

Table ICll gives the parameters obtained with this 
model. These results are used to compare the different mor- 
phologies assumed in the flat and flared models. 



